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We map the ground-state ensemble of antiferromagnetic 
Ising model of spin-5 1 on a triangular lattice to an interface 
model whose entropic fluctuations are proposed to be de- 
scribed by an effective Gaussian free energy, which enables 
us to calculate the critical exponents of various operators in 
terms of the stiffness constant of the interface. Monte Carlo 
simulations for the ground-state ensemble utilizing this inter- 
facial representation are performed to study both the dynami- 
cal and the static properties of the model. This method yields 
more accurate numerical results for the critical exponents. By 
varying the spin magnitude in the model, we find that the 
model exhibits three phases with a Kosterlitz-Thouless phase 
transition at | < Skt < 2 and a locking phase transition at 
| < Sl < 3. The phase diagram at finite temperatures is also 
discussed. 

PACS numbers: 75.I0.Hk, 64.60. Fr, 68.35. Rh, 75.40.Mg 



I. INTRODUCTION 

Over the years it has been found that there exist many 
two-dimensional classical spin models, discrete and con- 
tinuous alike, whose ground-state manifolds are macro- 
scopically degenerate and, more interestingly, also exhibit 
critical behaviours, i.e., spin-spin correlation functions 
within the ground-state ensembles decay with distance 
as power laws. The classification of universality class for 
these models has always been a challenging problem [jjj 

An earlier example of this kind is the antiferromag- 
netic Ising model on the triangular lattice. The exact 
solution for this model by Stephenson || showed that 
although this model remains paramagnetic at nonzero 
temperature, its ground state is critical. Later works 
by Blote et al revealed yet another remarkable property 
of the ground-state ensemble of this model, namely, it 
permits a Solid-on-Solid (SOS) representation in which 
spin fluctuations are subsequently described by the fluc- 
tuating interface in the SOS model g]. Recent stud- 
ies also demonstrated that this interfacial representation 
provides a valuable avenue for studying the ground-state 
ordering of quantum magnets |j|,D and the ground-state 
roughness of oriented elastic manifolds in random media 
Q. Other recently studied models with critical ground 
states include three-state antiferromatic Potts model on 
the Kagome lattice fjj,||, the 0(n) model on the hon- 
eycomb lattice P,[iT|, the Four-Coloring model on the 



square lattice [ |il]|l2| ] , and the square-lattice non-crossing 
dimer model and dimer-loop model |i~3[ ]. On the other 
hand, some very similar models with degenerate ground 
states exhibit long-range order, such as the constrained 
4-state Potts antiferromagnet . 

In this article we study the ground-state properties of 
antiferromagnetic Ising model of general spin on a tri- 
angular lattice which also belongs to the class of models 
mentioned above. Recent numerical studies of this model 
include Monte Carlo simulations Jl||,[l(| and transfer ma- 
trix calculations |r^| . Here we revisit this model by per- 
forming Monte Carlo simulations. The motivation of the 
present work is two-fold: (l)unlike previous simulations, 
we utilize the interfacial representation directly in ana- 
lyzing the simulation results, for example, we compute 
the stiffness constant of the fluctuating interface which, 
in turn, yields more accurate critical exponents of various 
operators; and (2) we also study the dynamical proper- 
ties of this model for the first time making use of the 
interfacial representation. 

The body of the this paper is organized as follows. 
Section [0| describes the model Hamiltonian and maps it 
onto a spin-1 problem whose interfacial representation is 
then described. In Section III, we propose an effective 



continuum theory for the long-wavelength fluctuations of 
the interface. Here we also show how to relate scaling 
dimensions of various operators to the stiffness constant 
of the interface, and derive some other analytical results 
based on this "height representation." This allows an- 
alytical understanding of the phase diagram (Sec. |rv| ). 
Details of Monte Carlo simulations and numerical results 
on both dynamical and static properties are presented in 
Section [v], including a comparison of the new and old ap- 
proaches to determining the exponents. As a conclusion, 
the paper is summarized and various possible extensions 
are outlined, in Section VI. 



II. THE MODEL 

The antiferromagnetic Ising model of spin-S* on a tri- 
angular lattice can be described by the following Hamil- 
tonian: 



(i) 



i 



where the spin variable s(r) defined on lattice site r of 
the triangular lattice can take any value from a discrete 
set [— 5, — 5+1, 5— 1, 5], and the sum over e runs over 
three nearest-neighbor vectors ei , e2 and e3 as shown in 
Fig. [I]. Here the coupling constant J is positive describing 
the antifcrromagnetic exchange interaction between two 
nearest-neighbor spins: s(r) and s(r + e). 

One important reason for interest in this model is that 
the 5 — > oo limit [p| is the same as the Ising limit of 
the (classical or quantum) Heisenberg antiferromagnet 
on the triangular lattice with Ising-like anisotropic ex- 
change. That model was shown to exhibit a continuous 
classical ground state degeneracy and unusual features of 
the selection by fluctuations of ground states Q . 

The ground-state configurations of the above model 
given by Eq. ([!]) consist of entirely of triangles on which 
one spin is +5, another is —5, and the third can be any- 
thing in [— 5, +5]. Thus, if spin s(r) takes an intermedi- 
ate value —S < s(r) < 5, this forces the six surrounding 
spins to alternate +5 and —5; exactly which interme- 
diate value s(r) takes does not matter in determining 
whether a configuration is allowed. 



A. Spin-1 mapping 

Therefore, this allows us to reduce each state {s(r)} to 
a state {cr(r)} of a spin-1 model, by mapping s(r) = +5 
to cr(r) = +1, s(r) = — 5 to cr(r) = — 1, and intermedi- 
ate values —5 < s(r) < +5 to er(r) = 0. In this spin-1 
representation of the model, the rules for allowed config- 
urations are exactly the same as for the 5 = 1 model; 
however instead of being equal, the statistical weights 
have a factor 25—1 for each spin with er(r) = 0. It 
should be noted that in the 5=1/2 case, s(r) = ±1/2 
simply maps to cr(r) = ±1. 

It can also be shown that the expectation of any poly- 
nomial in {s(r)}, in the ground-state ensemble of the 
spin-5 model, can be written in terms of a similar expec- 
tation in the spin-1 model. Specifically, one must simply 
replace 



s(rY 



S rn a{r), 

S m [(l-C m {S))a{vf 



m odd 
C m (5)] , m even 



(2) 

where (e.g.) C*2(5) = 1(1 — 5 _1 ). Thus there is no loss of 
information in this mapping. Indeed, in some sense, the 
extra freedom to have s(r) vary from —(5 — 1) to 5 — 1 
is trivial: once given that s(r) and s(r') are intermediate 
spin values, there is no correlation between these values. 

So we henceforth restrict ourselves to the spin-1 
mapped model whose partition function for its ground- 
state ensemble can be written as: 



Z= Y, (2S- 1 )' 



where n s denotes the number of zero spins in a ground- 
state configuration {er(r)}. By varying the weight factor 
continuously in the spin-1 model, it would possible to give 
a precise meaning to any real value of 5, and to simulate 
such an ensemble. However, in this article we perform 
Monte Carlo simulations for an ensemble in which 25 
takes only integer values. 

The spin-1 representation could be further reduced 
to a spin-1/2 representation cr(r) as described in 
Refs. |0|Jl3@. They let 



cr(r) = cr(r) + fc(r) 



(4) 



Here fc(r) = if a(r) = ±1 and if a(r) = 0, fc(r) = 
+1 or —1 according to whether the surrounding spins 
are (+1, — 1, +1, — 1, +1, — 1) or the reverse. Note this 
mapping is not invertible. The spin-1/2 representation is 
less satisfactory in that is arbitrarily breaks the up-down 
symmetry of correlation functions, but it was desirable 
for the transfer-matrix calculations of Lipowski et al J!?]] 
since it reduced the number of degrees of freedom. 



B. Height mapping 

We define a microscopic, discrete-valued height func- 
tion z(r) living on the vertex of the triangular lattice 
such that the step in z(r) between adjacent vertices is a 
function of the adjacent spins: 



z(r + e) — z(y) = — I — er(r + e)cr(r 



(5) 



(3) 
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where tr(r) is the spin-1 operator and e can be any of the 
three nearest-neighbor vectors e! 2,3- It is easy to show 
that the total change in height function, when traversed 
along any smallest loop, i.e, an elementary triangle, is 
zero. Therefore, z(r) is well-defined everywhere for the 
ground-state configurations, but it is not well-defined in 
any excited state. This prescription generalizes that orig- 
inally introduced by Blote et al for the case 5 = 1/2 
P, pT|j22| (the prescriptions agree in that case). 

This type of height mapping differs from other sorts 
of mapping (e.g. dualities) in a crucial way: since the 
spin microstates of the spin-1 model are mapped essen- 
tially one-to-one to the height microstates, it is possible 
to perform Monte Carlo simulations and construct con- 
figurations z(r) after each time step. We have found 
that analysis of the z(r) correlations is much more effi- 
cient for extracting critical exponents than analysis of the 
spin correlations directly as was done in previous Monte 
Carlo simulations |l~5| ]. 

III. HEIGHT REPRESENTATION THEORY 

In this section we propose an effective continuum the- 
ory which describes the long-wavelength fluctuations of 
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the interface. We also demonstrate how the critical expo- 
nents of various operators are determined by the stiffness 
constant of the interface. 



A. Effective free energy 

To describe the interface in the rough phase, we must 
define a smooth height field h(x) by coarse-graining the 
discrete field z(r). As a first stage, on every triangular 
plaquette formed by sites ri,r2,r3, define a new discrete 
height 



MR) = i(z(n)+z(r 2 ) + z(r 3 )) 



(6) 



where R is the center of a triangle. The possible values 
of the h(R) are {n/2}, for any integer n. (For the case 
S = 1/2, the only possible values are integers.) To each 
of these values corresponds a unique ground-state spin 
configuration of the spin-1 model on that triangle, i.e., 



s(r) = $ s (/i(r + u)-/i (r)) 



(7) 



where u is any vector from a site to the center of any 
adjoining triangle. The mapping is many-to-one: the 
function §> s (h) has period 6. Notice that the r.h.s. of 
Eq.(^) turns out to be independent of u, but the pe- 
riodic dependence on h is phase-shifted by a function 
ho(r) which takes different values on each of the three 
V3 x v^3 sublattices. Essentially, we have mapped the 
T = ensemble of the spin-1 problem into an equiva- 
lent interface problem. Note that, given a configuration 
of {/i(R)}, each a(r) is specified (via Eq. ([?])), once for 
each adjoining triangle. The requirement that these six 
values of cr(r) coincide translates to a somewhat compli- 
cated set of contraints between pairs /i(R) and /i(R') on 
adjoining triangles; the difference h(R) — h(R!) may be 0, 
±1/2, or ±1, but some of these are disallowed (depend- 
ing on which h() values are integer or half-odd-integer, 
and on the orientation of R — R'). The weight of each 
configuration is given, as in (||), by by (2S — 1)™ S . 

Fig. [I] shows the /i(R) mapping explicitly where the 
spins cr(r) take values from {+1,0,-1}. The twelve 
states are arranged in a circle because the pattern re- 
peats when h — > h ± 6. 

There are certain special "flat states" in which h(R) is 
uniform on all triangles. Each of these is periodic with 
a V3 x \/3 unit cell - in effect it is a repeat of one of 
the triangles in Fig. [p. We shall name these states by 
writing the spins on the three sublattices, "(+, +, — )" and 
"(+,-,0)"; here "±" stands for a = ±1. It should be 
noted that there are two non-equivalent species of flat 
state corresponding to integer, and half-integer valued 
/i(R) respectively. They are non-equivalent in the sense 
that they are not related by lattice symmetries. One of 
the species that is favored by the locking potential (see 




) below) is what is previously called "ideal" states 



Thus we can imagine that all states can be described as 
domains of uniform h(R) separated by domain walls. Fi- 
nally, by coarse-graining h(R) over distances large com- 
pared to the lattice constant, one obtains /i(x) which 
enters the conjectured continuum formula for the free 
energy, which is entropic in origin j^], 



F({fc(x)} = jd* 



f |vMx)| 2 



V(h(x)) 



(8) 



where K is the stiffness constant of the fluctuating inter- 
face. 

A lattice shift by one lattice constant leaves the free 
energy invariant, but induces global shifts in height space 
/i(x) — > /i(x) ± 1; hence the potential V(-) in (||) must 
have period one. It is typically approximated as 



V(h) « h v cos(2nh). 



(9) 



Such a periodic potential, usually referred as the locking 



term |23|, favors the heights to take their discrete values 
one of the two types of flat state, depending on the sign of 
hy- For large S we expect hy < 0, favoring the (+, — , 0) 
states, in view of the large entropy of flippable spins; it is 
not so sure which state is favored at smaller S. but this 



doe s not m atter for the critical exponents (see Sec. Ill C 
and III D| , below. 



B. Fluctuations and correlation functions 

In the rough phase, by definition, the locking term is ir- 
relevent, and so the long- wavelength fluctations of height 
variable /i(x) are governed by the Gaussian term of Eq. 



F({h(x)} = |dx||vMx)| 2 = ^|q 2 |Mq)l 2 , (io) 

where we have performed the Fourier transform. Hence 
by equipartition theorem, 



5,(q)^(|Mq)| 2 



1 



Aq 2 



(11) 



Similarly, we can also measure the height-height differ- 
ence function in the real space as: 

c,(R)^i([MR)-Mo)] 2 ) 

ln(7ri?/a) + ... (i?> 1) , (12) 



2ttA 



where a is the lattice spacing cutoff. 
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C. Scaling dimensions 



Using Eq. (12), we can compute the scaling dimension 
xo of any local operator 0(r), which is defined as in the 
correlation function, 



(O*(r)O(0))~f 



-1x 



(13) 



By local operator, we mean that O(r) is a local function 
of spin operators in the vicinity of r. Now, the same 
spin configuration is recovered when the height variable 
h(R) is increased by 6. |^4| Thus any local operator 0(r) 
is also a periodic function in the height space, and can 
consequently be expanded as a Fourier series: 



O(r) 



G 



G e 



iGh(r) 



JGqHt) 



(14) 



for spin and bond-energy operators respectively. 

If a magnetic field is implemented by adding a term 
— H ^ r <r(r) to the Hamiltonian, then our dimension- 
less uniform "magnetic field" is defined by H' = H/T. 
The exponents associated with H' (and with the uniform 
magnetic susceptibility), are easily related to the corre- 
lation exponents of the uniform magnetization operator, 



M(R) = i( < r(n) + <7(r 2 ) + < 7(r3)) 



(19) 



where R is the center of a triangle formed by sites 
r 1; r 2 , r 3 . A simple inspection of Fig. [[] shows that such 
an operator has a periodicity of 2 in the height space, 
thus yielding: 



27T 

(jM — "7T 



(20) 



where G runs over height-space reciprocal-lattice vectors 
(i.e. multiples of 27r/6). The last step of simplification 
in ( |l4| ) follows because the scaling dimension xo of the 
operator O(r) is determined by the leading relevant op- 
erator in the above expansion, i.e., Go is the smallest G 
with nonzero coefficient in the sum. Inserting Eq. Jl4| ) 
into Eq. (|l^) and making use of Eq. ([l2]) , we obtain the 
following: 



(0*(r)0(O)) = ( e -'GoMr) e .Gok(0)) 

= P - G o c h( r ) r -no 



(15) 



Therefore, the critical exponent r\o associated with the 
operator O(r) is given by: 



770 = 2ico = 



1 



2ttK 



\Go\ 2 



(16) 



E. Zone-corner singularities 

Observe that the microscopic height variable z(r) in 
any flat state is not uniform but is rapidly modulated 
with the wave vector Q = 4^(1,0). The amplitude of 
modulation itself is a periodic function of the coarse- 
grained height field /i(x) which in turn implies that the 
correlation function decays with distance as a power-law, 
and consequently that its structure factor has a power- 
law singularity at Q. 

Such a zone-corner singularity is also directly con- 
nected to the singularity in the structure factor of the 
bond-energy operator. To see this, recall that there is a 
linear relation between the microscopic height variables 
and the bond-energy operator given by Eqs. (^J) and jl^), 
i.e., 



D. Definition of operators 

In this paper, besides the usual spin operator <r(r), we 
also study the bond-energy operator E(r + e/2) for the 
reason that will become clear in the next section: 



E{r + e/2) = I + l a ( T + e)a(r) 



(17) 



where e denotes one of the three nearest-neighbor vectors 
as before. 

As discussed already, the spin operator on a given site 
has a periodicity of 6 in the height space, from which a 
simple inspection shows that the bond-energy operator 
is also periodic in the height space with a periodicity of 
3. Therefore, the reciprocal lattice vectors of the most 
relevant operator in the Fourier expansion in Eq. (H4) 



Ga 



2tt 

IT' 



G t 



2tt 



(18) 



E(r + — ) = z(r + e) — z(r) 



(21) 



Then it is interesting to note that the Fourier transform 
E e (q) of bond-energy operator given above turns out to 
be 

E e (q) = A" 1 / 2 V e l ^ r+ ^E{r + -) 



-2i sin(-q • e>(q) 



(22) 



In other words, as a byproduct of measuring (|z(q)| 2 ), 
we have at the same time measured the structure factor 
of, say, the bond-energy operator of the same orientation 
specified by the nearest-neighbor vector e: 

^(q)^(|£; e (q)| 2 )=4sin 2 (iq.e)(|z(q)| 2 > . (23) 



We will utilize this relation in Sec. |yPj to extract the 
exponent of bond-energy operator from the Monte Carlo 
simulations. 
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F. Exact solution for S = 1/2 

The 5 = 1/2 triangular Ising antiferromagnet is ex- 
actly soluble, by the same techniques which solve the 
ferromagnetic two-dimensional Ising model, and was im- 
mediately recognized to have critical behavior as T — > 0. 
The spin and energy correlation functions were computed 
exactly by Stephenson; it transpires that rj a — 1/2 and 
r]E = 2 exactly, implying through the arguments of Blote 
et al (see Sec. [II C and [II D| ) that the effective stiff- 
ness in Eq. (JTc ) is K = 7r/9 exactly. The exponents 
implied by the interface scenario || - in particular, the 
magnetic field exponent r/M - are fully confirmed by nu- 
merical transfer-matrix computations. [ p5| | 

The Coulomb gas picture of Kondev et al [^6| , wherein 
the 5=1/2 triangular Ising antiferromagnet is viewed 
as a fully-packed loop model with fugacity 1, also 
predicts the exact exponents. 



IV. PHASE DIAGRAM 

In this section, we collect some consequences of the 
height representation for the phase diagram and the na- 
ture of the various phases within it. |p8| 

A. Kosterlitz-Thouless and locking transitions 

The locking potential V(-) in (||) favors the flat states. 
In view of (||), its leading reciprocal-lattice vector is 
Gy — 27r, corresponding to a scaling index xy = 
\Gy\ 2 /irK = ir/K for the corresponding conjugate field 
hy. It is well known that if 2 — xy > 0, then hy becomes 
relevant (under renormalization and the interface locks 
into one of the flat states. |2^] Since K grows monotoni- 
cally with S, such a locking transition occurs at a critical 
S L where K L = tt/2 = 1.57079... |Jl7). In this "smooth" 
phase, any spin operator 0(r) has long-range order, by 
arguments as in Sec. |III C . 



have at most short-range (exponentially decaying) corre- 
lations. Thus we expect them to have a spatial "white 
noise" spectrum: 



B. Fluctuations in smooth phase 

One of our aims in this paper was to pinpoint the lock- 
ing transition Sl, which demands that we have a crite- 
rion to distinguish these phases. We must supplement 
Eq. (|ll|), which shows the expected qualitative behavior 
of height fluctuations (|ft.(q)| 2 ) in the rough phase, with 
a parallel understanding of the smooth phase. 

In the smooth state, the symmetry (of height shifts) is 
broken and a fully equilibrated system has long-range or- 
der, such that (/i(x)) is well defined and uniform through- 
out the system. Fluctuations around this height, then, 



(IMq)l 5 



const 



(24) 



for small q. 

A phase with "hidden" order was suggested by 
Lipowski and Horiguchi jl7],|2(|. Numerical transfer- 
matrix calculations ]l7j ] using the spin- 1/2 representation 
indicated < r\ a < 1/9 for 25* > 6, which is impossible 
if the spin correlations are derived from height fluctua- 
tions, |3fl as we reviewed in Sec. III. An exotic phase to 



reconcile these facts was to postulate a phase in which 
the interface was smooth and (5"(r)) ^ 0, yet for the real 
spins (c(r)) = as suggested by spin correlation mea- 
surements. 

What does this imply for our height variable h(R), 
which has a one-to-one correspondence with the real spin 
configuration {<r(r)}? If the interface is smooth, then 
the probability distribution of height values on a given 
plaquette, P(/i(R)), is well defined. In order to "hide" 
the order, it is necessary that P(h) correspond to zero 
expectations of the spins. Now, reversing s(r) on all three 
sites in the plaquette requires h — ► h ± 3, as seen from 
Fig. 0. One can convince oneself that, to have ensemble 
average (c(r)) = 0, the distribution P(h) must be at least 
as broad as ^5(h — hi) + ^5(h—fi2), with hi — h 2 = (h)±3, 
implying the bound 

Var[h(R)] = (h(Ti) 2 ) - (h(Ti)) 2 > (3/2) 2 . (25) 



C. Finite temperature behavior 

At T > 0, plaquettes with non-minimal energy are 
present and they correspond to vortices in the function 
/i(x). Thus, unfortunately, the height approach of ana- 
lyzing simulations more or less breaks down. Neverthe- 
less, one can still predict the T > phase diagram from 
knowledge of the T = stiffness constant derived from 
our simulations. The shape of this phase diagram has al- 
ready been explained in Ref. |jl7[| ; here we note some ad- 
ditional interesting behaviors which can be predicted (fol- 
lowing Ref. ||(b)) using the exponents associated with 
vortices. 

The other exponents in Kosterlitz-Thouless (KT) the- 
ory are associated with elementary defects (often called 
vortices). Indeed, it is easy to check (in this system) 
that the excess energy of a non-ground-state plaquette 
is directly proportional to its vortex charge (a Burgers 
vector in height-space) , so the effect of nonzero tempera- 
ture is simply to make the vortex fugacity nonzero. The 
vortex exponent is r\ v = so as usual the vortex fu- 

gacity becomes relevant and defects unbind, destroying 
the critical state, at the KT transition defined by a spin 
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exponent taking the critical value r\ a = 1/4. If r\ a > 1/4 
at zero temperature, i.e. if < Kkt = 2tt/9 — 0.69813..., 
then defects unbind as soon as T > 0. Thus a zero- 
temperature KT transition occurs at Skt defined by 
K = K KT . Jl| 

Ref. ji~7| did not, however, address the critical ex- 
ponents of the correlation length £(T) and the spe- 
cific heat C(T) as a function of temperature, which 
are also controlled by vortex exponents. Naively, if 
the energy cost creating one vortex is E c , and if the 
minimum excitation is a vortex pair, then one would 
expect the low-temperature specific heat to behave as 
C(T) ~ exp(-2i; c /T) and at S = 1/2 this is indeed true 
|p7j . However, the renormalization group ||] shows the 
singular specific heat behaves as 

f(T) ~ y{Tf/i^) (26) 

where y(T) = exp(—E c /T) is the vortex fugacity; conse- 
quently when rj v < 2, the true behavior is 

C{T) ~ exp(-2£i/T) (27) 

with Ei = 2E c /(4 — rj v ) < E c . (Physically, part of the 
excitation energy is cancelled by the large entropy due to 
the many places where the vortex pair could be placed.) 
This behavior has been observed in the 3-state Potts anti- 
ferromagnet on the Kagome lattice (?|] , and should occur 
in the present system for all 5 > 1/2. 

D. Finite magnetic field 

It is interesting to consider the effect of a nonzero mag- 
netic field H' . It is known already that at 5 — 1/2, ||] 
such a field is an irrelevant perturbation, so that the sys- 
tem remains in a critical state, yet at sufficiently large 
H it undergoes a locking into a smooth phase, [|25| ap- 
proximated by any of the three symmetry-equivalent flat 
states of type "(+,+,—)" with magnetization 5/3 

As also already noted jl7| , there is a critical value S c h 
defined by ^(Sch) =4/9, beyond which r\ m = 9t7 ct < 4 
so that the system locks into long-range order as soon as 
H' is turned on. Within this regime, there are still two 
subregimes with different behavior of M(h) near h = 0. 
For 2 < i]m < 4, the initial slope is zero, i.e., the sus- 
ceptibility is not divergent; when tjm < 2, as occurs for 
5 > 2, there is a divergent susceptibility and correspond- 
ingly there should be a singularity at q = in the spin 
structure factor (|er(q)| 2 ). 

What do we expect in the locked phase at S > 5l? 
Here the difference between the two kinds of flat states 
becomes crucial. The H' field favors the (+,+,—) type 
of flat state, but entropy favors the (+,—,0) type of 
flat state. Thus we expect a transition to the (+, +, — ) 
state only at a nonzero critical field H' c . On reducing 
H' through H' c , a twofold symmetry breaking occurs, in 



which one of the + sublattices becomes the (disordered) 
sublattice; hence, this transition should be in the Ising 
universality class. Presumably the line H' C (S) meets the 
H' = axis at S — Sl- There must also be line of lock- 
ing transitions S c h{H'), which terminates on the H' = 
axis at S c h- 

For 5 — 1/2, the effect of the magnetic field was con- 
firmed numerically in Ref. p{|. 

V. MONTE CARLO SIMULATIONS AND 
RESULTS 

In this section we describe the implementation details 
of Monte Carlo simulations performed for spin-1 model 
in which 25* takes only integer values from 1 to 8. We 
then present numerical results for the relaxation times 
of slow modes in the Monte Carlo dynamics. Two dif- 
ferent methods of compute the critical exponents of the 
spin, bond-energy, and uniform-magnetization operators 
are described in different sub-sections: one in terms of 
the extrapolated stiffness constants of the interface and 
the other in terms of the singularities of the corrsponding 
structure factors. 

A. Details of Monte Carlo Simulations 

A spin is called Qippable if its six surrounding nearest- 
neighbor spins alternate between +1 and —1. Clearly, 
changing the value of this flippable spin results in an- 
other new spin configuration in the ground-state ensem- 
ble, provided that we start with a spin configuration in 
the ensemble. Moreover, such an update maintains the 
global tilt of the interface due to the local nature of this 
update. This update will be used as our Monte Carlo up- 
date in this paper. Two slightly different cases arise for 
different values of 25*: (1) for 2S = 1, the local update is 
precisely equivalent to a spin flip i.e., er(r) — > — cr(r) due 
to the absence of zero spin; and (2) for all other values of 
2S, a random choice must be made in the local update: 
for example, cr(r) = — > cr(r) = 1 or —1. (Recall S 
denotes the spin magnitude of the original model.) 

Let n s and rif denote the number of zero-spins and 
flippable spins of configuration </>. If an attempted single- 
spin update for (f> results in a new configuration 4>' with n' s 
and n'j,, then the transition probability W in accordance 
with the detailed balance principle is: 

W= W -min{l,^-}-min{l,(2S-l) n °- n °} , (28) 

where Wo denotes the bare transition probability: Wo = 
J- for 25 = 1, and W = ^- for 25 > 2 which reflects 
the random choice to be made in the local update as dis- 
cussed above. With the transition probability given in 
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Eq. (p8|), it is straightforward to show that the detailed 
balance principle is satisfied, i.e., P((f))W((f> — > <j>') = 
P(<fr')W((f>' — > (f>), where P(4>) denotes the probability for 
configuration <fi to occur and P(4>) ~ (25— 1)™ S since each 
spin configuration in the original spin-5 model has equal 
probability to occur. Note also that nf/n'f = 1 + 0(1/N) 
for large N, so this rule is important only because of the 
finite system size. 

To implement in practice the transition probability 
given above, we randomly select a site out of a list of 
the n/ flippable sites, and randomly update this spin to 
one of the two possible new spin values if 25 > 2 or sim- 
ply flip this spin if 25 = 1. The total numbers of zero 
spins n' s and flippable spins in the resulting configura- 
tion are then computed. This update is subsequently ac- 
cepted with a probability: mm{l,nj/rij} • min{l, (25 — 

l)" s_ " s }. A practical implementation of the transition 
probability given in Eq. is thus achieved. 

Throughout this paper, a unit time or one Monte Carlo 
Sweep (MSC) is defined such that there are N s attempts 
of updating within this unit of time (or one attempt per 
spin on average). Here N s denotes the total number of 
spins in the simulation cell. The simulation cell always 
contains N s = 72 x 72 spins in this paper unless explicitly 
mentioned otherwise. Periodic boundary conditions are 
adopted. Since we always start with a flat state, the 
simulations are thus performed in the sector with a zero 
global tilt of the interface. 



B. Dynamical scaling: the relaxation time r q 

We now discuss the correlations between the configu- 
rations generated sequentially in the Monte Carlo simu- 
lations by studying the relaxation time of the slow modes 
in the model, namely, the Fourier modes /i q which play 
the role of an order parameter j^]. The linear-response 
dynamics of such a mode is usually formulated as a 
Langevin equation, 



d/i(x, t) 
dt 



Sh{x) 



■fat) 



(29) 



where T is the dissipation constant, and the static free 
energy functional F({/i(x)}) is given by Eq. (|J). Here 
£(x, t) is a stochastic noise generated in the Markov chain 
of Monte Carlo simulations. As it is expected that the 
correlation time of the slow mode under consideration is 
much longer than that of the noise, and since the update 
steps are local and independent, it is proper to model 
£(x, t) as Gaussian noise, uncorrelated in space or time: 



(fat) fa', t')) = 2T8(x-x')5(t-t') 



(30) 



in which the choice of 2r ensures that the steady-state 
of the interface under the Langevin equation (|2^) agrees 
with its equilibrium state under the free energy (B). 



This linear stochastic differential equation can be 
solved easily by performing Fourier transform. Eq. 
thus reduces to 



dh(q, t) 
Jt 



= -r^|q| 2 Mq,i) + £(q,i) 



(31) 



which implies an exponentially decaying correlation func- 
tion of (/i*(q, t)h(q, 0)) ~ e~ l l T ^ with the relaxation time 
r q given by 



TK l 



(32) 



Therefore, the dynamical scaling exponent for the Monte 
Carlo dynamics, defined by r q ~ |q| -2 , is always z = 2 
in the rough phase. 

To check this prediction on the dynamical scaling ex- 
ponent in practice where the above continuum theory is 
regularized on a lattice, we compute the following auto- 
correlation function C(q, i) of the microscopic height 
variable z(q): 



C(q,t) 



(z*(q,0)z(q,t))-|(z(q,0))| 2 
(z*(q,0)z( q) 0)>-Kz(q,0)>|2 



(33) 



Here () stands for the dynamical average, and the time 
t is measured in unit of MCS. For each interger-valued 
25 = 1, 2, 8, we perform 10 5 MCS's with a flat initial 
configuration and compute the auto-correlation functions 
upto t < 50 for modes that correspond to the five smallest 
|q| 2 values. In Fig. |^, we display the results so obtained 
for 25 = 1. Other cases of 25 are found to have very sim- 
ilar features. It is clear from Fig. |^ that log 10 C(q, t) can 
be fitted very well by a — i/r q where a and the relaxation 
time T q are the fitting parameters. In other words, the 
relaxation is strictly exponential in all cases. Note that 
we used a cutoff t = 10 in our fitting. The same fitting 
procedure is carried out for other cases of 25. 

The final results of the relaxation time T q as a func- 
tion of |q| 2 for 25 = 1, 6 are shown in Fig. 0; and for 
25 = 6, 7,8 as an insert. The fact that r q scales as |q| 2 
for 25 = 1, 5 as indicated by the fitting in Fig. || thus 
shows that the ground-state ensembles for 25 = 1, ...,5 
are in the rough phase. On the other hand, it is in- 
deed clear from the insert that for 25 = 7 and 8, r q 
curves downward as |q| 2 — > which is in sharp constrast 
to those of 25 — 1,...,5. From this, we conclude that 
ground-state ensembles for 25 = 7 and 8 are in the flat 
phase. As for 25 — 6, it is not conclusive from the data 
available whether r q scales as |q| 2 or curves downward 
as |q| 2 — *■ 0. Nonetheless, the fact that the relaxation 
time of the slowest mode for 25 = 6 is longer than for 
any smaller or larger value of 5, suggests that 25 = 6 is 
very close to the locking transition. Further support for 
this phase diagram is also obtained by explicit calcula- 
tions of stiffness constants and critical exponents which 
is discussed in the next section. 
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C. Stiffness constants and critical exponents 



zone corner q — > Q = ^(1,0) in the thermodynamic 
limit N s — > oo; and (2) for 25 — 1, it approaches a 



As implied by Sec. [II B , the stiffness constant of the 
fluctuating interface can be directly measured by study- 
ing the long-wavelength fluctuations of the height vari- 
ables, i.e., their structure factor as given by Eq. (fll|). It 
should be noted that concerning the task of calculating 
the Fourier components h(q) in Eq. (|TT|) , it can be re- 
placed by the approximation in terms of the microscopic 
height variables z(q) given by 



h(q) « z(q) 



WQ 



<z{r) 



(34) 



where r labels a lattice site of the finite triangular lat- 
tice of total N s lattice sites used in the simulation. Here 
wo = V3/2 is the weight of a lattice site, i.e., the area 
of its Voronoi region, which is introduced so that the mi- 
croscopic height variable z(q) coincides with the coarse- 
grained height variable h(q) in the long-wavelength limit 
(q — * 0). But unlike h(q), z(q) still contains fea tures 
such as zone-corner singularities discussed in Sec. Ill E 
that are only manifested in miscroscopic height variables. 

Starting with a flat state, we perform 2 x 10 3 MCS's as 
the equilibrium time; subsequent measurements of phys- 
ical quantities are carried out at intervals of 20 MCS's. 
This separation is a compromise between the correlation 
times of small q modes and of larger q modes, which are 
respectively much longer and somewhat shorter than 20 
MCS - see Fig. |. Each run consisted of 8 x 10 5 MCS, i.e. 
4 x 10 4 measurements were taken; these were subdivided 
into 20 independent groups for the purpose of estimat- 
ing statistical errors. The same procedure is used for all 
25 = 1,2,.. .,8 reported in this paper. 

In Fig. |, we plot (Iz(q)l 2 )- 1 vs. q 2 for 25 = 1, in- 
cluding all q in the first Brillouin zone. From the plot, 
we observe that (|^r(q) I 2 }" 1 is remarkably isotropic up to 
about q 2 ~ 1.5. This comes about because of the 6-fold 
rotational symmetry of the triangular lattice which en- 
sures that anisotropy occurs only in q e and higher order 
terms, assuming that the function is analytic. This is 
in constrast to other models defined on the square lat- 
tice where anisotropy already sets in at the order of g 4 
jij],[l0]]. The lower envelope of the data points in Fig. |] 
corresponds to the line of q y = in the g-vector space. 
Other cases of 25 are found to have very similar fea- 
tures as illustated in the insert of Fig. ^ where we plot 
the lower envelope for all 25 = 1, 2, 8. The structure 
factor of the height variables appears to diverge in the 
long- wavelength limit |q| 2 — > for all 5 values, even for 
the largest 5 values. (In the latter case, however, we 
believe one would see the plot asymptote to a constant 
value, in a sufficiently large system; see below.) 

Two other interesting features of the structure factor 
are also revealed in the insert in Fig. [|: (1) for 25 > 
2, it appears to indicate yet another singularity at the 



constant instead. As already discussed in Sec. HIE, the 



appearance of zone-corner singularities is expected, the 
precise nature of such singularities, however, is discussed 
in the next section. In the remaining of this section, 
we analyze the zone-center singularity to check if height 
variables behave as required by Eq. ([ll]) for the rough 
phase and consequently extract the stiffness constants. 

To further study the nature of zone-center singularity 
in terms of how (|z(q)| 2 ) scales as a function of q 2 in 
the long- wavelength limit, we show the log-log plot of 
(Iz(q)l 2 )- 1 vs. q 2 for 25 = 1,...,8 in Fig. § Compar- 
ing the simulation results for different systems sizes of 
L = 36, 48 and and 72, we notice that data are well con- 
verged down to accessible small q vectors - except for 
the case of 25 = 6 and 7, where the finite size effect is 
still discernible. This is, of course, consistent with the 
fact that 25 = 6 and 7 are close to the locking transition 
where the correlation length diverges; it is interesting, 
however, to notice that their finite-size trends are differ- 
ent. In the case 25 = 6, the data plot for L — 72 curves 
upwards less than that for L — 48, while in the case 
25 = 7, the L = 72 data show more upwards curvature 
than the L — 48 data. 

By fitting (|z(q)| 2 ) -1 to a function q 2a with a be- 
ing the fitting parameter, we obtain, using the data of 
system size L = 72 and a cutoff q 2 < 0.5, the expo- 
nent a = 0.990(1), 0.988(1), 0.986(2), 0.984(2), 0.974(2) 
and 0.935(1) respectively for 25 = 1,2,3,4,5, and 6. 
Apart from the case of 25 = 6, these values agree with 
a = 1 as in the predicted q -2 power-law singularity of 
the structure factor in the rough phase, Eq. ([11]) ■ As 
for 25 = 7 and 8, (|z(q)| 2 ) -1 clearly deviates from a 
power-law scaling and instead curves upwards to level 
off, which indicates that models with 25 = 7 and 8 are 
in the smooth phases wher e (|.z( q)| 2 ) remains finite as 



0, as discussed in Sec. IV B. This conclusion is in 



excellent agreement with that in ferre d from dynamical 
scaling analysis presented in Sec. VB . 

It should be noted that in Fig. |5|, as a general proce- 
dure adopted throughout this paper in extracting numer- 
ical values of some physical quantities, we have averaged 
the data corresponding to the same magnitude of |q| 2 to 
further reduce the effect due to statistical errors. The 
relative statistical error on each individual data point 
(|z(q)| 2 ) of small q, which is measured directly from the 
variance among the 20 groups, is found to range from 
1% to 3%. This is indeed consistent with the estimates 
of such relative errors from the relaxation times of the 
slowest modes of models with different values of 25 al- 



ready given in Sec. VB. It is perhaps also worth noting 
that another good check on the statistical errors on each 
data point is to compare the values of (|z(q)| 2 ) for three 
q vectors which are related by 120° rotations in recip- 
rocal space, which ought to be equal by symmetry. For 
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example, in the case of 25 = 1, the values of (|,z(q)| 2 ) 
for the three q vectors of the same smallest magnitude 
q 2 = 0.0101539 of system size L = 72 are, respectively, 
285.551, 280.528, and 280.566, from which one thus also 
obtain the relative error of about 1%. This observation 
therefore motivates the averaging procedure used in this 
paper. 

The stiffness constants can be subsequently determined 
by fitting q _2 (|z(q)| 2 ) -1 to the function K + Ciq 2 for 
the isotropic part of the data in which the stiffness con- 
stant K and C± are the fitting parameters. The final 
fitting on the averaged data is shown in Fig. ^| where we 
used a cutoff q 2 < 0.5 in the fitting. We also tried other 
different cutoffs of q 2 < 0.1 and q 2 < 1.0, and found as 
expected that the stiffness is not sensitive to the value 
of cutoff as long as it falls into the isotropic part of the 
data. For example, we obtain, in the case of 25 = 1, 
K = 0.3488±0.0022, 0.3490±0.0008, and 0.3488±0.0006 
for cutoff q 2 < 0.1,0.5, and 1.0 respectively. Therefore, 
taking into account of the uncertainty introduced due to 
the cutoff, our final estimate for the stiffness constant 
is then K — 0.349 ± 0.001 which is in excellent agree- 
ment with the exact value -Kexact = 0.349065.... Similar 
procedure is carried out for other cases of 25 and the re- 
sults are tabulated in Table I. In the same table, we also 
give the value for the critical exponents of spin, bond- 
energy and uniform magnetization operators which are 
obtained straightforwardly according to Eqs. ([l6|), fl8| ) 
and (20). The agreement of our r]& values with the 



"77,7" values from transfer-matrix eigenvalues (see Table I 
of Ref. J!?]] , is quite close and becomes better as 5 grows 
(until 25 = 6.) 

As already discussed in Sec. IV C, a Kosterlitz- 



Thouless (KT) transition occurs at a critical value Skt 
where r\ a = 1/4, such that for 5 > Skt algebraic cor- 
relations persist even at small finite temperatures. It is 
clear from our data that Skt > 3/2. 

As for 25 = 6, the value of q- 2 (|z(q)| 2 )- 1 = 1.75±0.06 
at the smallest nonzero q 2 = 0.010153 is already larger 
than K L = tt/2 = 1.57079. That is, even if the sys- 
tem may have a "rough" behavior at the length scales 
probed in the simulation, the stiffness constant is such 
that the locking potential is relevant and must domi- 
nates at sufficiently large length scales, as discussed in 
Sec. IV A. A similar observation has already been used 



to argue that the constrained Potts antiferromagnet is 
in a smooth phase jl4| . This fact together with the poor 
fitting using the formula suitable for the rough phase (see 
the top curve of Fig. ^) leave us little choice but to con- 
clude that the ground-state ensemble for 25 = 6 also 
falls into the smooth phase, or possibly is exactly at the 
locking transition. 

Just as the finite-size effect for 25 — 6 was severe both 
for the spin-spin correlations (measured via Monte Carlo 
J15|,^6|) and also in spin-operator eigenvalues (measured 



via tranfer-matrix, |jl7)) we similarly find it is severe for 
height fluctuations. However, in view of the exponential 
relationship between the exponents and the stiffness con- 
stant, the latter measurements are much more decisive as 
to the true phase of the system. 

To sum up, based on the analysis on the nature of 
the singularity in the height structure factor at the long- 
wavelength limit and the numerical results on the stiff- 
ness constants, we thus conclude that the model exhibits 
three phases with a KT phase transition at § < Skt < 2 
and a locking phase transition at | < 5l < 3. 



D. Structure factor and zone-corner singularity 

Another more traditional approach |l5| in calculating 
the critical exponents of various operators is to com- 
pute the corresponding structure factors and analyze the 
power-law singularities at the appropriate ordering wave 
vectors. Namely, if the correlation function of an opera- 
tor O decays with distance as power-law (thus critical) 



(O(r)O(O)) 



e iQr 

r"no 



(35) 



then its structure factor near the ordering vector Q shows 
a power-law singularity 



5o(q = Q + k) -k 2 ^ - 1 ) 



(36) 



from which the critical exponent r\o = 2xo can be numer- 
ically extracted. Here in this section, we adopt this ap- 
proach to calculate the critical exponents of spin, bond- 
energy, and uniform-magnetization operators so as to 
compare with those obtained from the stiffness constant. 

As given by Eq. (f|), 5 E (q=Q + k) ~ 
(|z(q = Q + k)| 2 ). Here Q =1^(1,0) is the ordering vec- 
tor of the bond-energy operator. Therefore the interest- 
ing feature of structure factor of height variables, namely, 
the appearance of zone-corner singularity as shown in 
Fig. ||, is not only expected but also very useful in ex- 
tracting the critical exponent t)e- 

Of course, such a zone-corner singularity can also be 
understood within the framework of interfacial represen- 
tation, as in Sec. III. particularly Subsec. HIE. (Similar 
zone-corner singularities have been studied in Refs. p"l| ] 
and |l3).) Finally, according to the exact result t]e = 2 
(xe = 1) in the case of 25 = 1, i.e., 5#(q = Q + k) ~ 
j c 2(a: E -i) _^ cons t ^ the puzzling absence of the zone- 
corner singularity for 25 = 1 as shown in Fig. || is also 
resolved. 

In Fig. 0, we plot log 10 5£;(q) vs. log 10 |q-Q| 2 
where we have averaged data points with the same 
magnitude of |q— Q| 2 . Fitting 5_e(q) to the function 
| q _ Q|2(**-i)(Ci+C 2 |q- Q|) where x E ,Ci and C 2 are 

the fitting parameters, we obtain the critical exponents 

(s) 

rj E which are tabulated in Table I. In practice, we used 
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two different cutoffs in the fitting: |q— Q| 2 < 0.1 and 
< 0.5. The fitting for the latter is shown in Fig. ^, and 
the final quoted errors take into account the uncertainty 
due to the cutoffs. 

Similarly, we also computed the structure factor for the 
spin operator 5<j(q) using Fast Fourier transform while 
computing the height-height correlation function within 
the same Monte Carlo simulations. Results are shown 
in Fig. U and the extracted exponents are also tabulated 
in Table I. Fitting precedure used is exactly the same as 
that for bond-energy except that we fit 5 CT (q) to the func- 
tion Ci|q — Q| 2 ( a " _1 ) with C\ and x a being the fitting- 
parameters. From Table I, we note that the critical expo- 
nents extracted in this way are in good agreement with 
those obtained from stiffness constant utilizing the inter- 
facial representation, however, the latter yields much bet- 
ter statistical errors by an order of magnitude using the 
same Monte Carlo simulation data. This clearly demon- 
states the superiority of the interfacial representation in 
extracting critical exponents from numerical data. Simi- 
lar points were made regarding other models, but based 
on much less extensive simulation data, in Refs. and 

IE!; 

Similar fits were attempted for 25 = 6, yielding 



6) = 0.53 ± 0.41 and <r)i S) {2S = 6) = 0.236 ± 



0.036. While the statistical error on 77^ (25 = 6) is too 
large to render the fitting meaningful, the increase in the 
value of i]i S \2S = 6) when compared with r]i S \2S = 5) 
is added evidence that 2S = 6 is not in the rough phase; 
if it were still rough at that value of 5, we would have 

(S) 

expected a continuation of the decreasing trend of iy a 
with 5. 

As for the cases of 25 = 7 and 8, the structure factors 
of both the spin and bond-energy operators show weaker 
than power-law behavior as q — ► Q, as in Figs. |^ and [|, 
but they increase to a larger value (not seen in these 
logarithmic plots) right at Q. This is indeed consistent 
with the 5-function singularity, expected if these cases 
fall into the smooth phase with long-ranged order of the 
spin and bond-energy operators. 

Finally, we consider the uniform magnetization cor- 
relation exponent t)m- When 5 > 3/2, it can be pre- 
dicted (see ry^ in Table I) that t)m < 2, implying a 
divergent (ferromagnetic) susceptibility and a divergent 
structure factor Sm (q) as q — > Now, due to the lin- 
ear relation fll9| ) between {Af(R)} and {er(r)}, we im- 
mediately obtain 5a/ (q) ~ 5<j(q) near q = 0, just 
as 5_e(q) ~ (|z(q)| 2 ) near q = Q (see Sec. HIE and 
Eq. (|2g|)) Thus, a singularity at q = is expected in 
the structure factor of spin operator which is plotted in 
Fig. |j. From this figure, it appears that only for 25 = 4, 
5, and 6 does Sm (q) show a power-law singularity indi- 
cated by a straight line in this log- log plot. This confirms 
the prediction based on the stiffness constant; however, 
the numerical values of t\m extracted this way (see Ta- 



ble I) differ considerably from those calculated from the 

stiffness constant in the case of 25 = 5 and 6. 

(s) 

It is also apparent from Table I that ry a is systemati- 
cally overestimated as compared with the more accurate 
value derived from height fluctuations. We suspect that a 
similar overestimation affected the values of r\ a that were 
deduced from the finite-size scaling of the susceptibility of 
the staggered magnetization |15|l(| (this obviously mea- 
sures the same fluctuations seen in 5er(q) near Q.) Those 
data (also quoted in Ref. ]T^]) have quoted errors about 
four times as large as ours for r)a . Their exponent val- 
ues are all noticeably larger than the accurate value (r^^ 
or r/oo from Ref. Jl7|) - becoming worse as 5 grows (for 
25 = 4, 5 the difference is twice their their quoted error.) 
Clearly the systematic contribution to their errors was 
underestimated. The transfer-matrix method |l7| ought 
to provide the effective exponent r\ a for spin correlations 
on length scales comparable to the strip width, and hence 
is likewise expected to overestimate r\ a ; indeed, every r\„ 
value found in Ref. is slightly larger than our corre- 
sponding j\a value. 



E. Smooth Phase 

Which type of flat state is actually selected in the 
smooth phase? Fig. [l(] shows the measured expectation 
of n Sl the number of zero spins in the spin-1 represen- 
tation, for 1 < 25 < 8. As 5 grows, it is found that 
(n s ) approaches its maximum allowed value N s /3 as in 
the (+, — , 0) state, rather than zero, as in the (+, +, — ) 
state. Thus, the flat states with half-integer valued h(TL) 
in Fig. [j] are being selected in the smooth phase. Trans- 
lating back to the spin-5 model, this means that spins on 
two sublattices of the triangular lattice take the extremal 
values, +5 and —5 respectively, while spins on the third 
sublattice remain disordered. 

It is perhaps more illuminating to study the distri- 
bution of height variables to probe the height fluctua- 
tions in the smooth phase. To this end, we also show, 
in Fig. |l0|, the histogram of height variable h(H) in the 
cases of 25 = 2 and 25 = 8, which is measured for a 
typical configuration generated in the Monte Carlo simu- 
lations. J3l[]. The broad distribution observed in the case 
of25 = 2(5<5i) evolves to a narrowly peaked dis- 
tribution in the case of 25 = 8 (5 > Sl)- (It decays as 
exp(— const|/i — (h)\).) Thi s supports the intuitive pic- 
ture presented in Sec. |VB. Furthermore, the center of 
this peaked distribution is half-integer valued. (Numeri- 
cally, the mean is (h) — 0.46 for the distribution plotted 
in Fig. [h].) In other words, the locking potential V(h) 
favors the (+, 0, — ) type of flat state, in which one sub- 
lattice is flippable, rather than the (+, +, — ) type of flat 
state. (See Fig. |l|). This kind of flat state was also ex- 
pected analytically in the limit of large 5 [p9 30 1. 
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We have also computed Var(/i) for each value of 5, 
in two ways. First, Var(z) is just normalization factors 
times Sq^od-^cOl 2 )' which we accumulated throughout 
the Monte Carlo run, as described earlier in this section; 
then it can be shown that Var(h) = Var(z) — 5 + \{n s ) 
exactly. For N s = 72 this gives Var(/i) = 1.06' and 0.20 
for 25 = 2 and 25 = 8, respectively, showing the contrast 
of the rough and smooth behavior. Secondly, we can 
compute Var(/i) directly from the histogram (from one 
snapshot) seen in Fig. |l^; this gives respective values 1.1 
and 0.15, in satisfactory agreement with the first method. 

The exotic "hidden order" phase ( see Sec - 



can be ruled out on the basis of these data: according to 
Eq. (25) the variance of /i(R) should be at least (3/2) 2 = 
2.25 in the hidden-order phase, while our measurements 
indicate it is at most only 0.20. Furthermore, for 25 = 
7 and 8, the structure factor S a (Q) at the zone-corner 
wave vector Q (not plotted) was much larger than at 
nearby q; that direct suggests a ^-function singularity in 
the thermodynamic limit, i.e., existence of long-ranged 
spin order in which (s(r)) ^ on at least two of the 
sublattices. 

Additionally, the spin structure factor S a (q) near the 
zone-corner wave vector Q (Fig. ||) showed a striking cur- 
vature in the "smooth" cases 25 = 7 and 8, quite differ- 
ent from the behavior at smaller 5. This makes it plau- 
sible that 5<j(q) — > constant, so that spin fluctuations 
have short-range rather than power-law correlations for 
5 > Sl- (It was not emphasized in Ref. p7| , but power- 
law correlations are implied if one takes seriously their 
measured values < 77^ < 1/9 for 25 = 7,8.) 

We propose, then, that actually rj a = i]e = i]m = 
for 5 > 5 C 2, as in the simplest picture of the smooth 
phase, and that the observed nonzero values are sim- 
ply finite-size effects due to the very slow crossover from 
rough to smoo th behavior near a roughening transition 
(see Sec. VI B , below, for further discussion.) 



VI. CONCLUSION AND DISCUSSION 

To conclude, in this article, we have investigated the 
ground-state properties of the antiferromagnetic Ising 
model of general spin on the triangular lattice by per- 
forming Monte Carlo simulations. Utilizing the inter- 
facial representation, we extrapolated the stiffness con- 
stants by studying the long-wavelength singularity in 
the height variables, which in turn lead to straightfor- 
ward calculation of critical exponents of various opera- 
tors within the framework of height representation. The 
results so obtained are further compared with those ex- 
tracted from a more tranditional method, and demon- 
strate that the method in terms of height representation 
method is by far the preferable one for extracting the 
critical exponents. We also analyzed both the dynamical 
and static properties of the model in order to map out 



the phase diagram which consists of three phases with 
a Kosterlitz-Thouless phase transition at | < Skt < 2 
and a locking phase transition at | < 5^ < 3. Even in 
the smooth state, analysis of the height fluctuations (as 
in Var(/i) was helpful in resolving questions which are 
made difficult by the strong finite-size effects near the 
locking transition. 



A. Rational exponents? 

One of our initial motivations for this study was the 
possibility of finding rational exponents even for 5 > 1/2. 
We believe the results in Table I are the first which are 
accurate enough to rule out this possibility. Indeed, 
7y CT (25 = 4) w 3/16 and n a (25 = 5) w 4/27, with dif- 
ferences similar to the error (0.001). But any random 
number differs from a rational number with denominator 
< 30 by the same typical error. The exception is that 
r]i K \2S = 6) is quite close to 1/9, but we have given 
other reasons to be suspicious of this value. 



B. What is 5 L ? 

Another intriguing question was whether the critical 
values 2Skt and 25l are exactly integers. Previous data 
Jl7t suggested that Sl = 3 exactly, and had large enough 
errors that Skt — 3/2 could not be excluded. Since 
Vo-(Skt) = 1/4 and tJo-^Sl) = 1/9, this question was 
answered by the preceding subsection: we find that def- 
initely Skt < 3/2. Furthermore, we suspect Sl < 3 as 



concluded in Sec. VC since the effective stiffness at the 



length scale we access is more than enough to drive the 
system to the locked phase. 

The question of the value of Sl suggests paying closer 
attention to the behavior of systems near the locking 
transition. It has been noted previously how the locked 
phase tends to behave qualitatively like the rough phase 
in a finite-size system, since the crossover is a very slow 
function of size. p5[] This is consistent with the appar- 
ent power-law behaviors observed at 5 > Sl in previous 
studies jl5|,[l7j and with the tendency of those studies to 
overestimate the exponents r]& and tje (as compared with 
our more accurate estimates.) This would suggest that, 
if extensive finite-size corrections were included in our 
analysis, they would reduce our estimate of Sl a bit fur- 
ther, i.e. we would more definitely conclude that 25 = 6 
is in the locked phase. 

Our analysis near the locking transition at Sl suffers 
from our ignorance of the expected functional form of the 
critical behavior as a function of 5 — Sl ■ A study of the 
roughening transition |32| used the Kosterlitz-Thouless 
(KT) renormalization group to derive analytic approx- 
imations for the total height fluctuation (closely analo- 
gous to Var(/i) in our problem), which made it possible to 
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overcome very strong finite-size effects and fit the rough- supported by NSF grant DMR-9214943. 
ening temperature precisely. Use of KT finite-size correc- 
tions was also essential in extracting meaningful numbers 
from transfer-matrix calculations near the locking tran- 
sition induced by a magnetic field in Ref. pj|. Thus, a 
similar adaptation of the KT renormalization group to 
give expressions for the behavior of (|z(q| 2 ), as a func- 
tion of (small) |q| and S — Sl, or the functional form 
of K(S) near Sl, could make possible a more conclusive 
answer as to whether Sl = 3 exactly. 



C. Possible improved algorithms 

Since the long-wavelength behavior in this model (in 
its rough phase) is purely Gaussian with z = 2 (see 
Sec. |V B| ) , the critical slowing down is particularly trans- 
parent. It seems feasible to take advantage of the exis- 
tence of a height representation to develop an acceler- 
ation algorithm. For example, it might be possible to 
extend the cluster algorithms which are known for the 
5 = 1/2 triangular Ising antiferromagnet. |33| These are 
well-defined at T > 0, but their effectiveness seems to 
depend in a hidden fashion on the existence of the height 
representation when T — > 0. 

An intriguing alternative approach starts from the ob- 
servation that at long wavelengths the system obeys 
Langevin dynamics (see Sec. VB and Ref. ||). Fourier 



acceleration |34j] , a nonlocal algorithm (efficient owing to 
use of the Fast Fourier Transform algorithm), is known to 
be effective in such cases. The key is to replace the uncor- 
related noise function £(x, t) of Eq. ( |30| ) with a new cor- 
related noise function having (|£(q, t)p) ~ 1 / 1 q| 2 . This 
might be implemented by first constructing a random 
function with such correlations, and then updating flip- 
pable spins with probabilities determined by that func- 
tion, in such a fashion as to satisfy detailed balance. 

Additionally, it may be possible to analyze transfer- 
matrices using the height representation. Quite possibly 
this would yield an order-of-magnitude improvement in 
the accuracy of the numerical results, for the same size 
system, similar to the improvement in analysis of Monte 
Carlo data. The transfer matrix breaks up into sectors 
corresponding to the step made by z(r) upon following a 
loop transverse to the strip (across the periodic boundary 
conditions. Then the stiffness can be extracted directly 
from the ratio of the dominant eigenvalues of two such 
sectors; such an analysis is already standard for qua- 
sicrystal random tilings, for which the long-wavelength 
degree of freedom is also an effective interface p5[. 
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TABLE I. Stiffness constant and critical exponents. Here rfc , V E an d Vm are ^ ne estimates for the critical exponents 
of spin and bond-energy operators calculated from the stiffness constant K as done in Sec. V(C), while r]i S \ rj^\ and Vm^ 
stand for the same critical exponents, but extracted from the singularities of their respective structure factors in Sec. V(D). 
Estimated errors are also given in the parenthesis. 



2S 


K 




Vh 


Vm 




vk 


(i') 

Vm 


1 


0.349(0.001) 


0.500(0.002) 


2.001(0.008) 


4.502(0.018) 


0.511(0.013) 


1.844(0.057) 




2 


0.554(0.003) 


0.315(0.001) 


1.260(0.006) 


2.836(0.013) 


0.332(0.016) 


1.340(0.072) 




3 


0.743(0.004) 


0.235(0.001) 


0.940(0.005) 


2.114(0.011) 


0.254(0.019) 


1.047(0.082) 




4 


0.941(0.006) 


0.186(0.001) 


0.742(0.004) 


1.670(0.010) 


0.203(0.022) 


0.791(0.092) 


1.634(0.014) 


5 


1.188(0.008) 


0.147(0.001) 


0.588(0.004) 


1.322(0.009) 


0.180(0.026) 


0.504(0.115) 


1.560(0.015) 


6 


1.597(0.015) 


0.109(0.001) 


0.437(0.004) 


0.984(0.009) 


0.236(0.036) 


0.530(0.410) 


1.527(0.016) 
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FIG. 1. Twelve flat states of the ground-state ensemble. 
Each flat state is simply specified by its spins on three sub- 
lattices A, B, and C of the triangular lattice since all spins 
on same sublattice take the same value. The height variable 
h(TVj defined at the center of an elementary triangle according 
to Eq. (4), which is uniform for each of these twelve states, is 
also shown. Note that /i(R) — > /i(R) + 6 results in identical 
spin configurations. The three nearest-neighbor vectors ei, 
e2 and e3 defined in Eq. (hi) are also displayed. 



FIG. 2. Auto-correlation function of Fourier components 
of the height variables -z(q) for 25* = 1. Only those corre- 
sponding to the five smallest |q| 2 are shown where we have 
averaged over data points with the same value of |q| 2 . The 
discrete q- vectors comes about because of the periodic bound- 
ary condition used for the Monte Carlo simulation cell which 
consists of 72 x 72 spins. The solid lines in the figure are the 
fittings discussed in Sec. IV(B) to extrapolate the relaxation 
time T q where we have used a cutoff in time t < 10 measured 
in the unit of MSG 



FIG. 3. Relaxation time r q as a function of q 



for 



25 = 1, 8. The solid lines are the fittings for cases of 
25 = 1, 6 from bottom to top (see Sec. V(B)). The dotted 
lines in the insert are only the guide for eyes where data for 
25 — 6, 7 and 8 are displayed from top to bottom. 



FIG. 4. Structure factor 5h(q) of height variables. We 
show in the main figure the inverse of the structure factor 
as a function of q 2 for 2S = 1. The lower envelope of the 
data corresponds to q y = 0. As an insert to the figure, we 
plot all the lower bounds for 2S = 1,2, ...,8 which go from 
bottom to top. Solid lines in the insert are only the guide for 
eyes. Note that q 2 = 17.5459... corresponds to the corner of 
the first Brillouin zone, i.e., q = Q = 4^(1,0). 



FIG. 5. Scaling of (^(q)! 2 )" 1 as a function of q 2 . We 
have averaged data points of the same magnitude of q-vector 
in each case of 2S = 1,2, ...,8 obtained for system sizes 
L = 36, 48, and 72. Note the error bars are smaller than 
the symbol size. Solid lines are fits using a cutoff q 2 < 0.5 
discussed in Sec. V(C). Dotted lines are only guides for the 
eye. 



FIG. 6. Extrapolation of stiffness constants. We show 
[q 2 (|z(q)| 2 )] _1 vs. q 2 as log-linear plot for 2S = 1,2,..., 6. 
Note that we have performed an average over data points with 
the same magnitude of q-vector for eac h ca se of 2S. Solid lines 
are the linear fitting discussed in Sec. V C in order to extract 



the stiffness constant which is given by the intercept of the 
fitting. Also note that the fittings shown are performed with 
a cutoff q 2 < 0.5. Fittings with other cutoffs are discussed in 
the text. 



FIG. 7. Structure factor Se(ci) of the bond-energy oper- 
ator near the zone corner Q. Data points are averaged re- 
sults over those with the same q — Q| 2 value for each case of 
25 = 1, 2, 8. Note that data points for each 25 have been 
shifted upwards by 0.5 with respect to their counterpart for 
25 — 1 in order to disentangle the data. Solid lines are the 
fittings discussed in Sec. V(D) to extract the critical expo- 
nent tje of the bond-energy operator. Dotted lines are only 
to guide the eye. 



FIG. 8. Structure factor S a (q) of the spin operator near the 
zone corner Q. Data points are averaged results over those 
with the same |q — Q| 2 value for each case of 25 = 1, 2, 8. 
Note that data points for each 25 are moved downwards by 
0.1 with respect to their counterpart for 25 — 1 in order to 
disentangle the data. Solid lines are the fittings discussed 
in Sec. V(D) to extract the critical exponent rj a of the spin 
operator. Dotted lines are only to guide the eye. 



FIG. 9. Structure factor 5a/ (q) of the spin operator near 
the zone center q — > 0. Data points are averaged results over 
those with the same |q| 2 value for each case of 25 = 1, 2, 8. 
Note that data points for each 25 are moved upwards by 
0.15 with respect to their counterpart for 25 — 1 in order to 
disentangle the data. Solid lines are the fittings discussed in 
Sec. V(D) to extract the critical exponent t]m of the uniform 
magnetization operator. Dotted lines are only to guide the 
eye. 



FIG. 10. Height distribution and ensemble average of the 
number of free spins n a , from one snapshot for each 5 value. 
On the top panel, we show the histograms of the height vari- 
ables h(R) for 25 = 2 and 25 = 8. On the bottom panel, 
n s is displayed as a function of 25. Note that the maximum 
allowed value for (n s ) is N s /3 where N s denotes the total 
number of spins in the simulation cell. 
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